%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                         %
%                          NeQuick-G - Test Code                          %
%                                                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%DESCRIPTION
%This short code is meant to show how to use the NeQuick-G function. The
%data were taken from the reference [1] test cases in Annex E.
%
%Reference: "Ionospheric Correction Algorithm for Galileo Single Frequency Users"
%Link: https://www.gsc-europa.eu/sites/default/files/sites/all/files/Galileo_Ionospheric_Model.pdf
%
%---------------------------- Pre-Settings --------------------------------

clear all; clc;

%Add Folders' Paths
addpath(genpath('./Geometr'));
addpath(genpath('./Ionosphere'));
addpath(genpath('./tables'));
addpath(genpath('./TEC'));

%-------------------------- Data Preparation ------------------------------

brdc_data = load("../../data/brde/Galieo_IONO.2019");
out_dir = "./data";
year = 2019;
nrows = 71;
ncols = 73;

gims = zeros(365*24, nrows, ncols);
for record = brdc_data'
    doy = record(1);
    [mo, dy] = ddd2mmdd(year, doy);
    %Time Data
    data.mth   = mo;  %Month [mth]
    %NeQuick-G Data Extraction of MODIP and F2/Fm3 Coeff.s Tables
    [data.stModip, data.F2, data.Fm3] = DataExtract(data.mth);
        
    %Broadcasted NeQuick-G Az Coeff.s (High Solar Activity)
    data.a0 = record(2); %236.831641;
    data.a1 = record(3);  % -0.39362878;
    data.a2 = record(4);   % 0.00402826613;
    %Galileo L1-Frequency [Hz]
    data.f      = 1575.42e+06;  
    for ut = 0:1:23
        idx_time = (doy-1)*24+ut+1;
        data.UT    = ut;  %Time  [UTC]
        ilat = 0;
        outfile = sprintf("%s/ne_%04d_%03d_%02d.grid", out_dir, year, doy, ut);
        grid = zeros(nrows,ncols);
        for lat = -87.5:2.5:87.5
            fprintf("doy:%03d, lat:%6.1f\n", doy, lat);
            ilat = ilat +1;
            ilon = 0;
            for lon = -180:5:180
                ilon = ilon +1;
                %lat = -45;
                %lon = -180;
                %Receiver Data
                data.RX_h   =  0;   %        78.11e-3;   %Altitude  [km]
                data.RX_lat =  lat; %          82.49;    %Latitude  [deg]
                data.RX_lon =  lon; %wrapTo180(297.66);   %Longitude [deg]
                %Transmitter Data
                data.TX_h   =  20000; %  19953735.00e-3;    %Altitude  [km]
                data.TX_lat =  lat;   %       35.85;     %Latitude  [deg]
                data.TX_lon =  lon;    % wrapTo180(140.68);    %Longitude [deg]
                
                %Test Case STEC [STECU]
                % STEC_real   = 37.81;
                
                %-------------------------- NeQuick-G Routine -----------------------------
                % [TEC] = TECVert(data.RX_h, data.TX_h, lat, lon, data);
                % disp(sprintf("RX: %.4f %.4f %.4f", data.RX_h, data.RX_lat, data.RX_lon));
                % disp(sprintf("TX: %.4f %.4f %.4f", data.TX_h, data.TX_lat, data.TX_lon));
                [TEC, ~, ~] = NeQuickG(data.RX_h, data.RX_lat, data.RX_lon,...
                                       data.TX_h, data.TX_lat, data.TX_lon,...
                                       data, data.f) ;
                grid(ilat, ilon) = TEC;
  
            end
        end
        save(outfile, "grid", "-mat");
        disp(outfile + " is saved!");
        gims(idx_time) = grid;
    end
end
save("./data/ne_2019.mat", "gims");
